home *** CD-ROM | disk | FTP | other *** search
/ CD Ware Multimedia 1995 May / cd Ware (Juegos) Epimundo.iso / DOS / C / CMATH.ZIP / MTST.C < prev    next >
Encoding:
C/C++ Source or Header  |  1988-12-08  |  9.1 KB  |  412 lines
  1. /*   mtst.c
  2.  Consistency tests for math functions.
  3.  Inverse function tests and Wronksians for random arguments.
  4.  With NTRIALS=10000, the following are typical results for
  5.  IEEE double precision arithmetic:
  6.  
  7. Consistency test of math functions.
  8. Max and rms relative errors for 10000 random arguments.
  9. x =   cbrt(   cube(x) ):  max = 1.54E-016   rms = 6.50E-018
  10. x =   atan(    tan(x) ):  max = 2.83E-016   rms = 7.68E-017
  11. x =    sin(   asin(x) ):  max = 4.44E-016   rms = 7.01E-017
  12. x =   sqrt( square(x) ):  max = 1.57E-016   rms = 2.70E-017
  13. x =    log(    exp(x) ):  max = 1.88E-014   rms = 2.02E-016
  14. x =   tanh(  atanh(x) ):  max = 4.00E-016   rms = 8.29E-017
  15. x =  asinh(   sinh(x) ):  max = 2.04E-016   rms = 1.56E-017
  16. x =  acosh(   cosh(x) ):  max = 2.67E-014   rms = 2.95E-016
  17. x = pow( pow(x,a),1/a ):  max = 3.38E-012   rms = 3.39E-014
  18. Legendre  ellpk,  ellpe:  max = 1.55E-011   rms = 1.74E-013
  19. x =  ndtri(   ndtr(x) ):  max = 4.09E-014   rms = 6.40E-016
  20. Absolute error criterion (but relative if >1):
  21. lgam(x) = log(gamma(x)):  max = 4.49E-016   rms = 8.80E-017
  22. x =    cos(   acos(x) ):  max = 4.44E-016   rms = 1.00E-016
  23. Absolute error and only 2000 trials:
  24. Wronksian of   Yn,   Jn:  max = 7.12E-011   rms = 1.59E-012
  25. Wronksian of   Kn,   Iv:  max = 7.07E-012   rms = 5.25E-013
  26.  
  27. */
  28.  
  29. /*
  30. Cephes Math Library Release 2.1:  December, 1988
  31. Copyright 1984, 1987, 1988 by Stephen L. Moshier
  32. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  33. */
  34.  
  35.  
  36. #include "mconf.h"
  37.  
  38. #define NTRIALS 1000
  39. #define WTRIALS (NTRIALS/5)
  40. #define STRTST 0
  41.  
  42. double ndtr(), ndtri(), ellpe(), ellpk(), gamma(), lgam();
  43. double fabs(), sqrt();
  44. double cbrt(), exp(), log(), tan(), atan();
  45. double sin(), asin(), cos(), acos(), pow();
  46. double tanh(), atanh(), sinh(), asinh(), cosh(), acosh();
  47. double jn(), yn(), iv(), kn();
  48.  
  49. /* Provide inverses for square root and cube root: */
  50. double square(x)
  51. double x;
  52. {
  53. return( x * x );
  54. }
  55.  
  56. double cube(x)
  57. double x;
  58. {
  59. return( x * x * x );
  60. }
  61.  
  62. /* lookup table for each function */
  63. struct fundef
  64.     {
  65.     char *nam1;        /* the function */
  66.     double (*name )();
  67.     char *nam2;        /* its inverse  */
  68.     double (*inv )();
  69.     int nargs;        /* number of function arguments */
  70.     int tstyp;        /* type code of the function */
  71.     long ctrl;        /* relative error flag */
  72.     double arg1w;        /* width of domain for 1st arg */
  73.     double arg1l;        /* lower bound domain 1st arg */
  74.     long arg1f;        /* flags, e.g. integer arg */
  75.     double arg2w;        /* same info for args 2, 3, 4 */
  76.     double arg2l;
  77.     long arg2f;
  78. /*
  79.     double arg3w;
  80.     double arg3l;
  81.     long arg3f;
  82.     double arg4w;
  83.     double arg4l;
  84.     long arg4f;
  85. */
  86.     };
  87.  
  88.  
  89. /* fundef.ctrl bits: */
  90. #define RELERR 1
  91. #define EXPSCAL 4
  92.  
  93. /* fundef.tstyp  test types: */
  94. #define POWER 1 
  95. #define ELLIP 2 
  96. #define GAMMA 3
  97. #define WRONK1 4
  98. #define WRONK2 5
  99. #define WRONK3 6
  100.  
  101. /* fundef.argNf  argument flag bits: */
  102. #define INT 2
  103.  
  104. extern double MINLOG;
  105. extern double MAXLOG;
  106. extern double PI;
  107. extern double PIO2;
  108. /*
  109. define MINLOG -170.0
  110. define MAXLOG +170.0
  111. define PI 3.14159265358979323846
  112. define PIO2 1.570796326794896619
  113. */
  114.  
  115. #define NTESTS 15
  116. struct fundef defs[NTESTS] = {
  117. {"  cube",   cube,   "  cbrt",   cbrt, 1, 0, 1, 2000.0, -1000.0,   0,
  118. 0.0, 0.0, 0},
  119. {"   tan",    tan,   "  atan",   atan, 1, 0, 1,    0.0,     0.0,  0,
  120. 0.0, 0.0, 0},
  121. {"  asin",   asin,   "   sin",    sin, 1, 0, 1,   2.0,      -1.0,  0,
  122. 0.0, 0.0, 0},
  123. {"square", square,   "  sqrt",   sqrt, 1, 0, 1, 170.0,    -85.0, EXPSCAL,
  124. 0.0, 0.0, 0},
  125. {"   exp",    exp,   "   log",    log, 1, 0, 1, 340.0,    -170.0,  0,
  126. 0.0, 0.0, 0},
  127. {" atanh",  atanh,   "  tanh",   tanh, 1, 0, 1,    2.0,    -1.0,  0,
  128. 0.0, 0.0, 0},
  129. {"  sinh",   sinh,   " asinh",  asinh, 1, 0, 1, 340.0,   0.0,  0,
  130. 0.0, 0.0, 0},
  131. {"  cosh",   cosh,   " acosh",  acosh, 1, 0, 1, 340.0,      0.0,  0,
  132. 0.0, 0.0, 0},
  133. {"pow",       pow,      "pow",    pow, 2, POWER, 1, 25.0, 0.0,   0,
  134. 50.0, -25.0, 0},
  135. {" ellpe",  ellpe,   " ellpk",  ellpk, 1, ELLIP, 1, 1.0,    0.0,  0,
  136. 0.0, 0.0, 0},
  137. {"  ndtr",   ndtr,   " ndtri",  ndtri, 1, 0, 1,   10.0,   -10.0,   0,
  138. 0.0, 0.0, 0},
  139. /* Absolute error criterion starts with tstyp = GAMMA: */
  140. {"gamma",  gamma,   "lgam",   lgam, 1, GAMMA, 0, 11.0, 1.0,   0,
  141. 0.0, 0.0, 0},
  142. {"  acos",   acos,   "   cos",    cos, 1, 0, 0,   2.0,      -1.0,  0,
  143. 0.0, 0.0, 0},
  144. /* Smaller number of trials starts with tstyp = WRONK1: */
  145. {"  Jn",     jn,   "  Yn",     yn, 2, WRONK1, 0, 30.0,  0.1,  0,
  146. 40.0, -20.0, INT},
  147. {"  Iv",     iv,   "  Kn",     kn, 2, WRONK2, 0,  30.0, 0.1,  0,
  148. 20.0, 0.0, INT},
  149. };
  150.  
  151. static char *headrs[] = {
  152. "x = %s( %s(x) ): ",
  153. "x = %s( %s(x,a),1/a ): ",    /* power */
  154. "Legendre %s, %s: ",        /* ellip */
  155. "%s(x) = log(%s(x)): ",        /* gamma */
  156. "Wronksian of %s, %s: ",
  157. "Wronksian of %s, %s: ",
  158. "Wronksian of %s, %s: "
  159. };
  160.  
  161. static double y1 = 0.0;
  162. static double y2 = 0.0;
  163. static double y3 = 0.0;
  164. static double y4 = 0.0;
  165. static double a = 0.0;
  166. static double b = 0.0;
  167. static double c = 0.0;
  168. static double x = 0.0;
  169. static double y = 0.0;
  170. static double z = 0.0;
  171. static double e = 0.0;
  172. static double max = 0.0;
  173. static double rmsa = 0.0;
  174. static double rms = 0.0;
  175. static double ave = 0.0;
  176.  
  177.  
  178. main()
  179. {
  180. double (*fun )();
  181. double (*ifun )();
  182. struct fundef *d;
  183. int i, j, k, itst;
  184. int m, ntr;
  185.  
  186. ntr = NTRIALS;
  187. printf( "Consistency test of math functions.\n" );
  188. printf( "Max and rms relative errors for %d random arguments.\n",
  189.     ntr );
  190.  
  191. /* Initialize machine dependent parameters: */
  192. defs[1].arg1w = PI;
  193. defs[1].arg1l = -PI/2.0;
  194. defs[3].arg1w = MAXLOG;
  195. defs[3].arg1l = -MAXLOG/2.0;
  196. defs[4].arg1w = 2.0*MAXLOG;
  197. defs[4].arg1l = -MAXLOG;
  198. defs[6].arg1w = 2.0*MAXLOG;
  199. defs[6].arg1l = -MAXLOG;
  200. defs[7].arg1w = MAXLOG;
  201. defs[7].arg1l = 0.0;
  202.  
  203.  
  204. /* Outer loop, on the test number: */
  205.  
  206. for( itst=STRTST; itst<NTESTS; itst++ )
  207. {
  208. d = &defs[itst];
  209. m = 0;
  210. max = 0.0;
  211. rmsa = 0.0;
  212. ave = 0.0;
  213. fun = d->name;
  214. ifun = d->inv;
  215.  
  216. /* Absolute error criterion starts with gamma function
  217.  * (put all such at end of table)
  218.  */
  219. if( d->tstyp == GAMMA )
  220.     printf( "Absolute error criterion (but relative if >1):\n" );
  221.  
  222. /* Smaller number of trials for Wronksians
  223.  * (put them at end of list)
  224.  */
  225. if( d->tstyp == WRONK1 )
  226.     {
  227.     ntr = WTRIALS;
  228.     printf( "Absolute error and only %d trials:\n", ntr );
  229.     }
  230.  
  231. printf( headrs[d->tstyp], d->nam2, d->nam1 );
  232.  
  233. for( i=0; i<ntr; i++ )
  234. {
  235. m++;
  236.  
  237. /* make random number(s) in desired range(s) */
  238. switch( d->nargs )
  239. {
  240.  
  241. default:
  242. goto illegn;
  243.     
  244. case 2:
  245. drand( &a );
  246. a = d->arg2w *  ( a - 1.0 )  +  d->arg2l;
  247. if( d->arg2f & EXPSCAL )
  248.     {
  249.     a = exp(a);
  250.     drand( &y2 );
  251.     a -= 1.0e-13 * a * y2;
  252.     }
  253. if( d->arg2f & INT )
  254.     {
  255.     k = a + 0.25;
  256.     a = k;
  257.     }
  258.  
  259. case 1:
  260. drand( &x );
  261. x = d->arg1w *  ( x - 1.0 )  +  d->arg1l;
  262. if( d->arg1f & EXPSCAL )
  263.     {
  264.     x = exp(x);
  265.     drand( &a );
  266.     x += 1.0e-13 * x * a;
  267.     }
  268. }
  269.  
  270.  
  271. /* compute function under test */
  272. switch( d->nargs )
  273.     {
  274.     case 1:
  275.     switch( d->tstyp )
  276.         {
  277.         case ELLIP:
  278.         y1 = ( *(fun) )(x);
  279.         y2 = ( *(fun) )(1.0-x);
  280.         y3 = ( *(ifun) )(x);
  281.         y4 = ( *(ifun) )(1.0-x);
  282.         break;
  283.  
  284.         case GAMMA:
  285.         y = lgam(x);
  286.         x = log( gamma(x) );
  287.         break;
  288.         
  289.         default:
  290.         z = ( *(fun) )(x);
  291.         y = ( *(ifun) )(z);
  292.         }
  293.     break;
  294.     
  295.     case 2:
  296.     if( d->arg2f & INT )
  297.         {
  298.         switch( d->tstyp )
  299.             {
  300.             case WRONK1:
  301.             y1 = (*fun)( k, x ); /* jn */
  302.             y2 = (*fun)( k+1, x );
  303.             y3 = (*ifun)( k, x ); /* yn */
  304.             y4 = (*ifun)( k+1, x );    
  305.             break;
  306.  
  307.             case WRONK2:
  308.             y1 = (*fun)( a, x ); /* iv */
  309.             y2 = (*fun)( a+1.0, x );
  310.             y3 = (*ifun)( k, x ); /* kn */    
  311.             y4 = (*ifun)( k+1, x );    
  312.             break;
  313.  
  314.             default:
  315.             z = (*fun)( k, x );
  316.             y = (*ifun)( k, z );
  317.             }
  318.         }
  319.     else
  320.         {
  321.         if( d->tstyp == POWER )
  322.             {
  323.             z = (*fun)( x, a );
  324.             y = (*ifun)( z, 1.0/a );
  325.             }
  326.         else
  327.             {
  328.             z = (*fun)( a, x );
  329.             y = (*ifun)( a, z );
  330.             }
  331.         }
  332.     break;
  333.  
  334.  
  335.     default:
  336. illegn:
  337.     printf( "Illegal nargs= %d", d->nargs );
  338.     exit(1);
  339.     }    
  340.  
  341. switch( d->tstyp )
  342.     {
  343.     case WRONK1:
  344.     e = (y2*y3 - y1*y4) - 2.0/(PI*x); /* Jn, Yn */
  345.     break;
  346.  
  347.     case WRONK2:
  348.     e = (y2*y3 + y1*y4) - 1.0/x; /* In, Kn */
  349.     break;
  350.     
  351.     case ELLIP:
  352.     e = (y1-y3)*y4 + y3*y2 - PIO2;
  353.     break;
  354.  
  355.     default:
  356.     e = y - x;
  357.     break;
  358.     }
  359.  
  360. if( d->ctrl & RELERR )
  361.     e /= x;
  362. else
  363.     {
  364.     if( fabs(x) > 1.0 )
  365.         e /= x;
  366.     }
  367.  
  368. ave += e;
  369. /* absolute value of error */
  370. if( e < 0 )
  371.     e = -e;
  372.  
  373. /* peak detect the error */
  374. if( e > max )
  375.     {
  376.     max = e;
  377.  
  378.     if( e > 1.0e-10 )
  379.         {
  380.         printf("x %.6E z %.6E y %.6E max %.4E\n",
  381.          x, z, y, max);
  382.         if( d->tstyp >= WRONK1 )
  383.             {
  384.         printf( "y1 %.4E y2 %.4E y3 %.4E y4 %.4E k %d x %.4E\n",
  385.          y1, y2, y3, y4, k, x );
  386.             }
  387.         }
  388.  
  389. /*
  390.     printf("%.8E %.8E %.4E %6ld \n", x, y, max, n);
  391.     printf("%d %.8E %.8E %.4E %6ld \n", k, x, y, max, n);
  392.     printf("%.6E %.6E %.6E %.4E %6ld \n", a, x, y, max, n);
  393.     printf("%.6E %.6E %.6E %.6E %.4E %6ld \n", a, b, x, y, max, n);
  394.     printf("%.4E %.4E %.4E %.4E %.4E %.4E %6ld \n",
  395.         a, b, c, x, y, max, n);
  396. */
  397.     }
  398.  
  399. /* accumulate rms error    */
  400. e *= 1.0e16;    /* adjust range */
  401. rmsa += e * e;    /* accumulate the square of the error */
  402. endlup:
  403.     ;
  404. }
  405.  
  406. /* report after NTRIALS trials */
  407. rms = 1.0e-16 * sqrt( rmsa/m );
  408. printf(" max = %.2E   rms = %.2E\n", max, rms );
  409. } /* loop on itst */
  410.  
  411. }
  412.